Unified theoretical framework for black carbon mixing state allows greater accuracy of climate effect estimation

Black carbon (BC) plays an important role in the climate system because of its strong warming effect, yet the magnitude of this effect is highly uncertain owing to the complex mixing state of aerosols. Here we build a unified theoretical framework to describe BC’s mixing states, linking dynamic processes to BC coating thickness distribution, and show its self-similarity for sites in diverse environments. The size distribution of BC-containing particles is found to follow a universal law and is independent of BC core size. A new mixing state module is established based on this finding and successfully applied in global and regional models, which increases the accuracy of aerosol climate effect estimations. Our theoretical framework links observations with model simulations in both mixing state description and light absorption quantification.

assuming either an internal or external mixture of aerosols 2,27 , leading to a wide range of estimated BC mass absorption cross-sections (MAC, a typical indicator of BC light absorption ability) from 3.1 to 18.0 m 2 /g (at 550 nm) on global average 27,28 . Therefore, a precise description of BC mixing state becomes the determinant factor of model performance when estimating BC optical properties and radiative forcing.
In this study, we built a theoretical framework linking dynamic processes to BC coating thickness distribution and discovered the selfsimilarity of BC mixing states, which was verified in eight different observation sites globally. The size distribution of BC-containing particles is found to follow a universal law. This self-similarity allows us to greatly simplify the characterization of BC mixing states in both model simulations and field observations. A new mixing state scheme was established for model simulation, which can precisely represent the BC mixing state. Model estimated BC absorption and radiative forcing is substantially reduced, which fits well with available field observation results. Our study links observations with model simulations in both mixing state and light absorption.

A universal law of BC mixing state
To characterize the BC mixing state, which is controlled by complex atmospheric processes, we performed a theoretical derivation considering the main physical processes affecting BC in the atmosphere. We discovered that the size distribution of BC-containing particles follows a universal law and is independent of BC core size. Figure 1 provides a conceptional scheme describing the main physical processes and evolution of aerosols in the atmosphere. We assume a monodisperse aerosol population (consisting of BC cores only) emitted into the atmosphere at time zero with diameter D c and number concentration n(D c ). After being emitted into the atmosphere, the particles experience both growth and deposition progresses continuously, which form a steady state 29 , that is, the size distribution of BC-containing particles is approximately steady (although the overall mass concentration may change). The growth of BC-containing particles via condensation and coagulation results in an enlarged particle size, whose change as a function of time is represented by the growth rate (GR).
The time evolution of the diameter of BC-containing particles (D p ) can be integrated to give At the same time, BC particles are removed by deposition process with the rate of Dep. The number concentration of particles at D p , i.e., n(D p ), decays due to deposition process is Then, the time evolution of n(D p ) can be integrated as Based on the steady-state approximation, Eqs. 2 and 4 can be combined and time term t can be eliminated. Taking the logarithm on both sides, we obtain the following equation: The slope k = Dep GR is determined by the deposition rate and the growth rate, and the intercept is determined by n(D c ). Note that we adopted a simplified derivation in the above theoretical analysis for better understanding. A more rigorous theoretical derivation as well as the interpretation of the dependency of GR and Dep on D p and time can be found in the Supplementary Information (SI). Equation 5 demonstrates that for different particle sizes, ln(n(D p )) and D p − D c (ΔD p , defined as coating thickness) are in a linear relationship (Fig. 1). The average coating thickness can be derived as 1/k (the detailed formula is shown as Eq. 7 in the "Methods" section). Furthermore, we find that the slope k is independent of D c , indicating the self-similarity of BC size distributions, that is, BC-containing aerosol with different core sizes should have similar distributions of coating thickness. Such self-similarity allows us to greatly simplify the description of BC mixing states since BC coating thickness can be fully described with known slope k and the BC core number-size distribution.
We verified our theoretical model by field observations of BC size using the single particle soot photometer (SP2, Droplet Measurement Technologies, USA), which is a well-recognized instrument to measure BC mixing state) from eight sites covering different environments globally 24,30,31 . As presented in Fig. 2, the BC size distribution follows an exponential law at all sites despite their different regions and properties (e.g., urban, regional background), which agrees well with our theoretical model. The slope k of linear regression, ranging from 0.008 to 0.020, provides a useful parameter to quantify BC size distribution and its absorption enhancement.
The D p distributions with different BC core sizes from SP2 measurements at four sites are presented in Fig. 3. We selected four D c bins (110-120 nm, 120-130 nm, 130-140 nm, and 140-150 nm) and calculated the average D p distribution in each bin. The result shows that the BC size distributions with different D c have approximately the same slope of ln(n(D p ))~D p . This phenomenon can be observed at all four sites, validating our theoretical model presented as Eq. 5, i.e., the shape of the BC size distribution is independent of D c . The presence of the same pattern in Nanjing (suburban), Maqu (remote background), Tokyo (urban), and at the Amazon Tall Tower Observatory site (affected by biomass burning) further indicates that the self-similarity of BC size distribution is ubiquitous in the real atmosphere.

Improved estimation of BC absorption and radiative effect
When evaluating BC absorption and radiative effect, simply assuming either internal or external mixture of aerosols may induce large discrepancies and uncertainty. Also, it is difficult to provide an accurate description of mixing state due to its great complexity. Based on our theory, the description of BC mixing state can be greatly simplified, which is applicable to the light absorption calculation in both climate models and chemical transport models. The absorption enhancement factor, E abs , is the ratio between the aerosol absorption coefficients before and after removal of coating, which is a widely used parameter to represent BC light absorption amplification. As shown in Fig. S3, the relationship between E abs and ΔD p is approximately linear when ΔD p is small. Therefore, the BC coating thickness distribution can be replaced by a monodisperse coating thickness, 1/k, when calculating the black carbon absorption (a detailed demonstration can be found in the "Methods" section). Light absorption coefficients over all D p and with mean D p show good agreement (Fig. S4), which further validates this simplification. Hence, this approximation can be applied directly in global and regional models for optical estimation.
Based on the above findings, we established a new mixing state module and applied it in a global climate model (CESM-CAM6) and a chemical transport model (WRF-chem) as examples. Model simulations of E abs and BC direct radiative forcing (DRF BC ) were performed using these two models alternatively with our module and the conventional assumption of mixing state (Figs. 4 and 5). Comparing the simulated E abs with observational data shown in Fig. 4, the simulated E abs using the conventional assumption of mixing state in both CESM-CAM6 and WRF-Chem (blue squares) are significantly higher than observations (1.0 to 1.7, shown as black squares). CESM-CAM6 uses the volume mixing assumption for E abs calculation and the simulated result is higher than 2.5. WRF-Chem has two types of mixing state assumption, which are volume mixing and core-shell mixing. However, the simulated E abs using both types of assumptions ranges from 2.0 to 2.5, which is also nearly twice that of observational results. Using the new module, the calculated E abs in CESM-CAM6 and WRF-Chem are 1.4 (1.3-1.6) and 1.4 (1.3-1.7), respectively, which agrees well with observational data, demonstrating the good performance of our mixing state description in the quantification of E abs .  Figure 5 presents simulated DRF BC using the new scheme and the conventional scheme in CESM-CAM6. We selected four typical regions (Europe, North America, South America, and Asia) to perform statistical analyses. Calculated DRF BC using our scheme is greatly reduced in all four regions (40-50%) compared with adopting the conventional scheme. This result is in agreement with the recently-found overestimation of radiative warming by BC in global climate models, largely due to the treatment of aerosol mixing state 32 . Figure 5 demonstrates that the new scheme based on our theoretical framework can efficiently resolve the existing overestimation of model-simulated radiative warming by BC.

Discussion
We built a unified theoretical framework based on a steady-state assumption to describe the complex mixing state of BC and find its self-similarity confirmed across a wide range of field observation sites. This universal law is the result of a balance of multiple atmospheric processes. Our findings link the representation of particle diameter (from field observations) and dynamic parameters (generally used in models), making observational data applicable in model simulations. Moreover, this unified theoretical framework reduces the dimension of mixing state descriptions. The mixing state module established in this study can be embedded in various types of atmospheric models and efficiently improves the accuracy of aerosol climate effect estimations without increasing computational complexity. We find that BC absorption enhancement and warming effect are much lower than current estimates.

Light absorption enhancement
We preformed the optical calculation using the core-shell Mie method. We used a lognormal distribution of D c . The geometric standard deviation was set to 1.8. The mean diameter of D c was 70 nm and the wavelength was 550 nm. The refractive indices (RI) of the BC and scattering components were set to 1.85 + 0.7i according to Bond et al. 1 and 1.53 + 0i according to Pitchford et al. 2 . To derive the response of MAC to ΔD p (Fig. S3), ΔD p from 10 nm to 200 nm with 10 nm intervals was adopted to calculate the mass absorption coefficient.
For the optical calculations in Fig. S4, we used the integral method (ΔD p varied from 1 nm to 1000 nm with 1 nm interval) and the k-value method to calculate absorption coefficients. The D p size distribution represented as Eq. 5 was adopted and k was set to 0.014. The setting of RI, wavelength, and D c size distribution were same as in Fig. S3.
As shown in Fig. S3, the relationship between E abs and ΔD p is approximately linear when ΔD p is small. Therefore, E abs could be represented as Eq. 6. where α(D c ) is the linear coefficient of E abs and ΔD p . With the known calculation formula of D p , we can derive that If Eq. 6 is integrated over all D p , the average E abs with given D c is found to be where c abs-external represents the light absorption coefficients of BC core. Equations 7 and 8 show that 1/k plays a similar role with coating thickness.

Field observations and site descriptions
Observational data of BC mixing states was collected from different sites, including Nanjing (a regional background site in the Yangtze River Delta region in China), Beijing (an urban site in the capital of China), Shaoguan (a regional background site in the Pearl River Delta region in China), the Tibetan Plateau (including three sites with different features), Japan (Tokyo), and the United States (Sacramento, influenced by biomass burning).
Field measurements in Nanjing were performed at the Station for Observing Regional Processes of the Earth system (SORPES, 118°57′10″ E, 32°07′14′′ N; 40 m a.s.l.), a regional background station in the western YRD region. A detailed description of SORPES can be found in previous studies 33,34 . Due to its geographical position, this observation platform is ideal to capture the transport from megacities in the YRD region and North China Plain. Observational data from February 2020, April 2020, and December 2021 was used in this study.
Observations in Lulang and Maqu were made over the Tibetan Plateau (TP) from April to July 2021. The Lulang site is located on the southeastern part of the TP with few traffic emissions nearby. The measurement period in Lulang was from 1 April to 25 May 2021. Maqu can be considered as a background site over the TP. The measurement

Conventional model
Modified model Observation Fig. 4 | Model simulated black carbon (BC) absorption enhancement using the new scheme in this study (red squares) and the conventional scheme (blue squares), compared with observations. WRF-Chem v and WRF-Chem cs stand for WRF-Chem simulations with volume mixing and core-shell mode, respectively. The black squares with error bars (standard deviation) are E abs observed using the thermodenuder (TD) method at different sites obtained from previous studies 9,13,[48][49][50][51][52] . There are two exceptions, which are Knox et al. 48 and Ueda et al. 50 . The error bars in Knox et al. 48 represent the E abs of aerosol with different age categories and the error bars in Ueda et al. 50 cover the 25th−75th percentiles. The E abs reported in Cappa et al. 51 were observed in two cities (Fresno and Fontana, California, respectively).  Table S1. Besides field observations conducted in this study, observational data from several sites covering different kinds of environment globally was also collected to support our findings. The measurement periods in Shaoguan, Beijing 31 , Tokyo 30 , Sacramento 9 , and the Amazon Tall Tower Observatory (ATTO) were December 2020, November 2014, August 2012, June 2010, and October 2019, respectively. Shaoguan is a regional background site in southeastern China and the Beijing site represents a heavily polluted region in the North China Plain. The Tokyo and Sacramento sites are located in urban areas in Japan and the United States, respectively. The selected observation period at the ATTO site covers a biomass burning episode including some relatively clean days. Detailed information on these observations is summarized in Table S1.
The physical properties of refractory BC particles were measured using single particle soot photometers (SP2, Droplet Measurement Technologies, USA). The operation principle of the SP2 has been well described in previous literature 20,35 . Briefly, sampled particles pass through a 1064 nm Nd:YAG laser beam. BC-containing particles heat up to their boiling point and incandesce. The BC mass can be computed based on its proportional relationship with the peak intensity of the incandescence signal and the BC mass equivalent diameter can be calculated with the known density of BC (normally assumed to be 1.8 g cm −3 36 ). The scattering calibration was performed using polystyrene latex spheres (PSL). The rBC mass was calibrated using fullerene soot with known diameter, which was selected by a differential mobility analyzer (DMA, TSI Inc., USA) and its mass was calculated using effective density values presented by Gysel et al. 37 . The leadingedge-only (LEO) fit method developed by Gao et al. 20 was adopted to calculate the scattering cross-sections of BC-containing particles and saturated scattering particles. Therefore, the optical diameters of BCcontaining particles can be further determined based on core-shell Mie theory.

Model simulations
There are two typical methods that are extensively applied for optical calculations in global climate models and regional transport models, i.e., the volume-mixing and core-shell Mie methods. The volumemixing Mie algorithm assumes that all components are mixed in all individual particles and the mean refractive index is calculated as the volume-weighted average of the refractive indices of each species. The core-shell Mie method assumes that BC is in the center, and other components are coated on the BC core. The shell refractive indices are assumed to be the volume-weighted average of the refractive indices of dissolved components. The volume-mixing Mie algorithm is included in CESM-CAM6. Both volume-mixing and core-shell Mie methods can be used to estimate aerosol optical properties in WRF-Chem. The refractive indices for shortwave radiation and densities of aerosol species in CESM-CAM6 and WRF-Chem model are summarized in Table S3. E abs is the ratio of MAC internal /MAC external , where MAC stands for the mass absorption cross section of BC. Since there is no external mixing module in CESM-CAM6 and WRF-Chem, the estimation of MAC external is based on off-line calculation. The default refractive indices and densities in each model are adopted. The BC diameter is assumed to follow a lognormal distribution with a count median diameter of 70 nm 2 .
The new mixing state and optical scheme based on our theory frame is established and applied in both CESM-CAM6 and WRF-Chem, which cover different model types. CESM-CAM6 and WRF-Chem are examples of global climate models and regional transport models, respectively. Moreover, CESM-CAM6 uses a modal aerosol module and WRF-Chem uses a sectional aerosol module, which are the two most widely implemented modules in atmospheric models. In the new scheme, a monodisperse coating thickness of 70 nm derived from k = 0.014 is adopted. In CESM-CAM6 and WRF-Chem, the BC core diameter is assumed to follow a lognormal distribution with a count median diameter of 70 nm.
We used the Community Atmosphere Model version 2.1.3 of the Community Earth System Model version 6 (CESM2.1.3-CAM6) 38 in the simulation of light absorption by BC and the global DRF with Modal Aerosol Module 4 (MAM4) 39 . MAM4 includes six aerosol components (BC, sea salt, sulfate, POA, SOA, and dust), which are divided into four modes (primary carbon mode, Aitken mode, accumulation mode and coarse mode), and simulates the mass mixing ratios of six components within each mode. The spatial resolution in the global simulation is 1.9°× 2.5°for a latitude and longitude grid with 70 vertical layers (from 50 m to~140 km). The simulation is performed for four years (2012-2015) with the spin-up in the first three years and analysis in the last year. The radiative transfer module in the shortwave is calculated by the radiation code RRTMG. The diagnostic calculation of CESM-CAM6 is conducted for the radiative properties of one specific component, namely by double running cases with and without that component. The aerosols in the accumulation mode in this study are resolved with a sectional representation (30 size bins) in the optical calculations.
The BC-induced direct radiative forcing (DRF BC ) in the conventional models is simulated using the default setting. The DRF BC determined with our new module by using k is performed assuming DRF BC is linear with MAC 2,40 . Thus, the DRF BC_k can be estimated from the change of MAC k and MAC in CESM.
WRF-Chem version 3.7 (Weather Research and Forecasting model coupled with Chemistry) was employed in this study, which is an online-coupled meteorology and chemistry model considering multiple physical and chemical processes, including emission and deposition of pollutants, advection and diffusion, gaseous 41 and aqueous chemical transformations, aerosol chemistry, and dynamics 42 . The model has been incorporated in several studies concerning the estimation of aerosol optical properties and its radiative forcing 43,44 . The model domain is centered at 35.0°N, 110.0°E with a grid resolution of 20 km that covers eastern China and the surrounding regions. A total of 30 vertical layers extending from the surface to 50 hPa are utilized in the model. The simulation is conducted for the first two weeks of April 2020, each run covers 36 h, and the last 24 h of output were kept for further analysis. The initial and boundary conditions of meteorological fields were updated from the 6-h NCEP (National Centers for Environment Prediction) global final analysis (FNL) data with 1°× 1°spatial resolution. The Rapid Radiative Transfer Model shortwave and longwave radiation scheme (RRTMG) represents the radiation transfer within the atmosphere 45 . Anthropogenic emissions from power plants, residential combustion, industrial processes, on-road mobile sources, and agricultural activities were derived from the MIX Asian emission inventory database 46 . Emissions of major pollutants such as carbon monoxide, sulfur dioxide, nitrogen oxides, ammonia, and speciated VOCs are all included. The MEGAN (Model of Emissions of Gases and Aerosols from Nature, version2) model embedded in WRF-Chem is used to calculate biogenic emissions online. Soil derived dust emissions are characterized by the GOCART emission schemes. For numerical representation of atmospheric chemistry, we used the CBMZ (Carbon-Bond Mechanism version Z) photochemical mechanism combined with the MOSAIC (Model for Simulating Aerosol Interactions and Chemistry) aerosol model 47 . Major aerosol components include BC, organic mass, sulfate, nitrate, ammonium, and other inorganic species. All aerosols were assumed to be spherical particles. The size distribution was divided into four discrete size bins defined by their lower and upper dry particle diameters (0.039-0.156, 0.156-0.625, 0.625-2.5, 2.5-10.0 μm). In each bin, aerosols were assumed to be internally mixed.

Data availability
The observation and simulation data generated in this study have been deposited in the Figshare database [https://doi.org/10.6084/m9. figshare.22490959]. Data collected from published papers are mentioned in the main text or the SI with corresponding references. Additional data related to this paper may be requested from the authors.

Code availability
The CESM2 source code can be downloaded from the CESM official website: http://www2.cesm.ucar.edu. The WRF-Chem source code can be downloaded from the WRF-Chem official website: https://www2. mmm.ucar.edu/wrf/users/download/get_source.html. Code related to this paper is available at https://doi.org/10.6084/m9.figshare.22455511. Additional code may be requested from the authors.